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Experimental measurements of the specific heat in glass-forming systems are obtained from the 
linear response to either slow cooling (or heating) or to oscillatory perturbations with a given fre- 
quency about a constant temperature. The latter method gives rise to a complex specific heat with 
the constraint that the zero frequency limit of the real part should be identified with thermodynamic 
measurements. Such measurements reveal anomalies in the temperature dependence of the specific 
heat, including the so called "specific heat peak" in the vicinity of the glass transition. The aim 
of this paper is to provide theoretical explanations of these anomalies in general and a quantitative 
theory in the case of a simple model of glass-formation. We first present new simulation results for 
the specific heat in a classical model of a binary mixture glass-former. We show that in addition 
to the formerly observed specific heat peak there is a second peak at lower temperatures which was 
not observable in earlier simulations. Second, we present a general relation between the specific 
heat, a caloric quantity, and the bulk modulus of the material, a mechanical quantity, and thus 
offer a smooth connection between the liquid and amorphous solid states. The central result of this 
paper is a connection between the micro-melting of clusters in the system and the appearance of 
specific heat peaks; we explain the appearance of two peaks by the micro-melting of two types of 
clusters. We relate the two peaks to changes in the bulk and shear moduli. We propose that the 
phenomenon of glass-formation is accompanied by a fast change in the bulk and the shear moduli, 
but these fast changes occur in different ranges of the temperature. Lastly, we demonstrate how 
to construct a theory of the frequency dependent complex specific heat, expected from heteroge- 
neous clustering in the liquid state of glass formers. A specific example is provided in the context 
of our model for the dynamics of glycerol. We show that the frequency dependence is determined 
by the same a-relaxation mechanism that operates when measuring the viscosity or the dielectric 
relaxation spectrum. The theoretical frequency dependent specific heat agrees well with experimen- 
tal measurements on glycerol. We conclude the paper by stating that there is nothing universal 
about the temperature dependence of the specific heat in glass formers - unfortunately one needs to 
understand each case by itself. 



I. INTRODUCTION 



The traditional measurements of the specific heat Cy 
at constant volume or Cp at constant pressure involve 
cooling (or heating) the sample at a constant rate [1, 2]. 
When applied to glass-forming systems, this approach 
has an inherent difficulty. Since glass-forming system 
tend to relax to equilibrium slower and slower as the tem- 
perature is lowered, at some point the 'constant rate' of 
cooling becomes too high for the system to respond to, 
and then the system does not reach equilibrium. Typi- 
cally the specific heat then drops abruptly, giving rise to 
the "specific heat peak" at some temperature which is 
sometimes identified as the glass transition temperature 
T g . Needless to say, such a definition of transition tem- 
perature is less than compelling, since it clearly depends 
on the rate of cooling and is not inherent to the system 
properties. 

In an attempt to overcome this difficulty Birge and 
Nagel [3] introduced " specific heat spectroscopy" . In this 
technique one keeps the sample close to a temperature T 
at all times, but perturbs it periodically with a small- 
amplitude oscillation of frequency lu. Linear response 
theory then relates the amount of heat exchanged at that 
frequency, SQ(lo) to the oscillatory temperature pertur- 



bation 5T(lu) via the relation 

6Q(u) = C{u)ST(u) 



(1) 



where C(u>) is the frequency dependent specific heat that 
can be measured at either constant volume or constant 
pressure. In order to find the thermodynamic specific 
heat one needs to extrapolate data to the w — > limit. 
Whether or not this extrapolation overcomes the above- 
mentioned worry of sufficient relaxation time is an issue 
that has not been fully clarified in the literature. 

In this paper we concern ourselves with the theoretical 
calculation of the specific heat in glass-forming systems 
and in the relation of the specific heat to other mate- 
rial properties. To this aim we focus on one simulational 
example (a binary mixture of point particles interacting 
via an r -12 repulsive potential) and one experimental 
example (glycerol). In the context of the first example 
we present results of new simulations that exhibit two 
distinct peaks in the curve of the specific heat vs. the 
temperature. We present for this example various the- 
oretical results, culminating with a new scenario to ex- 
plain the specific heat peaks, i.e. the micro-melting of 
clusters. We believe that this is the central point of the 
present paper. To understand the nature of the specific 
heat anomalies one must understand the physics that is 
behind the glassy behavior of this model in general and 
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FIG. 1: (Color online). A snapshot of the system at T = 0.44. 
In colours we highlight the clusters of large particles in local 
hexagonal order. The colours have no meaning. 



the existence of the two specific heat peaks in particu- 
lar. When the temperature is lowered at a fixed pressure 
this system [4] (as well as many other glass-formers [5-10] 
tends to form micro-clusters of local order. In the present 
case large particles form long-lived patches of hexagonal 
ordering first (starting at about T = 0.5, and at lower 
temperatures (around T — 0.1) also the small particles 
form long-lived hexagonal clusters. The clusters are not 
that huge, with at most 0(100) particles, (cf. Fig. 1), 
depending on the temperature and the aging time. But 
we have shown that the long time properties of correla- 
tion functions are entirely carried by the micro-clusters 
[4]. Below we will refer to the micro-clusters as curds 
and the liquid phase as whey. We will argue that the 
specific heat responds to the micro-melting of the clus- 
ters - those of small particles at the lowest temperatures 
and those of the larger particles at higher temperatures. 
The large increase in the number of degrees of freedom 
when a particle leaves a crystalline cluster and joins the 
liquid background is the basic reason for the increase in 
entropy that is seen as a specific heat peak. 

In the context of the second example we show that 
the calculation of the frequency dependent specific heat 
is easy when we have a reasonable model of the glassy 
relaxation. Having such a model for glycerol [11], we 
demonstrate in section IV that the information gained 
from the frequency dependent specific heat is very simi- 
lar to that learned from other linear response functions 
like broad-band dielectric spectroscopy. We will be able 
therefore to present spectra of the frequency dependent 
specific heat in close correspondence with experiments. 

The specific heat has interesting relations with the me- 



chanical moduli of the material, and we present relations 
(which pertain to any system with an r~" potential) to 
the bulk and shear moduli. As a result of our think- 
ing we conclude that the bulk and shear moduli change 
rapidly in the temperature range of the two distinct spe- 
cific heat peaks mentioned above. The relation to the 
bulk modulus is explicit, and is shown rigorously in Sub- 
sect. IIC. The relation to the shear modulus is less ex- 
plicit, simply because one does not have an equation of 
state with strain (a quantity that is ill defined in the 
context of glasses) . Nevertheless we present a conjecture 
that the bulk and shear moduli in generic glasses may 
change rapidly at two different temperature ranges. 

The structure of the paper is as follows: In Sect. II 
we present the model glass consisting of a binary mix- 
ture of point particles interacting via soft potentials, and 
discuss its thermodynamic properties. We derive exact 
equations for its specific heat at constant volume, which 
are correct at all temperatures through the glass transi- 
tion. We compare these results to molecular dynamics 
simulations in which special care had been taken to equi- 
librate the system, summarizing a computational effort 
of about two years. The main conclusion of this section 
is that the details of the interaction potential are crucial 
in determining what the specific heat does in the vicin- 
ity of the glass transition, and there is nothing universal 
about it. For a simple enough potential we can derive 
a theory that is in excellent agreement with simulations 
up to the first specific heat peak. To explain both peaks 
we must present a theory that takes into account explic- 
itly the tendency of the system to form micro-clusters 
[4]. The state of the system then becomes like curds of 
local crystalline order embedded in a whey of disordered 
fluid. It is the freezing or melting of these curds that 
account quantitatively for the specific heat peaks, as is 
shown in Subsect. III. Below we use interchangeably the 
words 'clusters' and 'curds'. In Sect. IV we turn to dis- 
cussing the frequency-dependent complex specific heat. 
To construct a theory of this object one needs a model 
of the dynamics of the system under study, be it glyc- 
erol or any other material. We demonstrate, using our 
dynamical model of glycerol [11], how this measurement 
is equivalent in terms of its dynamical contents to any 
other linear response to an oscillatory perturbation. We 
present theoretical spectra and show satisfactory agree- 
ment with the experiments. The paper ends in Sect. V 
where we draw conclusions and summarize the results 
and the implications of our calculations. 



II. THE BINARY MODEL AND ITS SPECIFIC 
HEAT 

The model discussed here is the classical example [12, 
14] of a glass-forming binary mixture of N particles in 
a 2-dimensional domain of area V, interacting via a soft 
1/r 12 repulsion with a 'diameter' ratio of 1.4. We refer 
the reader to the extensive work done on this system 



3 



[12, 14-17]. The sum-up of this work is that the model 
is a bona fide glass-forming liquid meeting all the criteria 
of a glass transition. 

In short, the system consists of an equimolar mixture of 
two types of particles, "large" with 'diameter' 02 = 1.4 
and "small" with 'diameter' <j\ = 1, respectively, but 
with the same mass m. In general, the three pairwise 
additive interactions are given by the purely repulsive — 
soft-core potentials 

Mr)=e(^fY , a, b =1,2, (2) 

where a aa = a a and a a b = {o~ a + crfc)/2. The cutoff radii 
of the interaction are set at 4.5tr a 6. The units of mass, 
length, time and temperature are m, ax, r = aiy/m/e 
and e/fcs, respectively, with ks being Boltzmann's con- 
stant. In numerical calculations the stiffness parameter 
of the potential (2) was chosen to be n = 12. 

We turn now to the analysis of the specific heat of this 
model as a function of the temperature. 

A. Specific heat (simulations) 




U/N 



FIG. 2: Color online: the energy distribution functions for the 
binary mixture model, computed at constant volume, such 
that the volume agrees with the pressure P=13.5 [14] at each 
temperature. 



The specific heat capacity (specific heat) at constant 
volume is defined by: 



Cv 
N 



d_{U)_ 
dT N 



(3) 



where d is the space-dimension and the potential energy 
of a binary mixture is given by: 



U 



1 



2^ 

in 



(4) 



Here is the distance between particles i and j. The 
average value of the potential energy is defined by aver- 
aging over configurational space T: 



(U) = 



J[/exp(-f )dr 
Jexp(-f )dr 



(5) 



Substitution of (5) into (3) yields the following expression 
for the specific heat: 



Cv 
N 



d (U 2 ) ~ {U? 
2 NT 2 



(0) 



The specific heat of our binary mixture model was mea- 
sured at constant volume in [13, 14, 18] and by us. In 
simulations one can measure the specific heat directly 
from its definition (3) or (6). We have used the last equa- 
tion which allows one to estimate the specific heat in a 
single run of the canonical ensemble Monte Carlo simu- 
lations. At each temperature the density was chosen in 
accordance with the simulation results in an NPT ensem- 
ble as described in [14] with the pressure value fixed at 



P = 13.5. As the initial configuration in the Monte Carlo 
process the last configuration of the molecular dynamics 
run for this model at given temperature after 1.3 x 10 8 
time steps was used. After short equilibration the po- 
tential energy distribution functions were measured dur- 
ing 2 x 10 6 Monte Carlo sweeps. The acceptance rate 
was chosen to be 30%. Simulations were performed with 
N = 1024 particles in a square cell with periodic bound- 
ery conditions. 

Examples of the spline interpolation of the potential 
energy distribution for a few temperatures are shown 
in Fig. 2. The first and second moments of these distri- 
butions define the average value of the potential energy 
(Fig. 3) and the specific heat (Fig4). We stress that these 
results were computed at constant volume, such that the 
volume corresponds to simulations in NPT ensemble with 
the pressure P=13.5 [14] at each temperature. 

One can see from these figures that the behavior of 
both quantities, the first and second moments of the dis- 
tribution, change abruptly in the vicinity of T ~ 0.5. The 
specific heat displays a maximum in the temperature de- 
pendence. Our simulations appear to provide trustable 
values of CV down to lowest temperatures where the 
value of the specific heat coincides with that of two- 
dimensional solid, i.e. Cv = 2. What could not be seen 
in earlier simulations is that there is a much smaller sec- 
ond peak of the specific heat at lower temperatures. To 
resolve it to the naked eye we present in Fig. 4 a blow- 
up of the region of lowest temperatures where the second 
peak is more obvious. To understand the nature of the 
specific heat anomalies we turn now to a theoretical anal- 
ysis of the caloric equation of state in order to study the 
specific heat using the definition (3). The physical origin 
of the two peaks will be explained in Subsect. III. The 
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FIG. 3: Color online: in dots: the temperature dependence of 
the average potential energy per in the binary mixture model, 
computed at constant volume, such that the volume agrees 
with the pressure P=13.5 [14] at each temperature. The con- 
tinuous line represents the approximation furnished by the 
virial expansion, which obviously fails for T < 0.5. 
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FIG. 4: Color online: in dots: the temperature dependence 
of the specific heat in the binary mixture model, computed at 
constant volume, such that the volume agrees with the pres- 
sure P=13.5 [14] at each temperature. The data indicate the 
existence of two specific heat peaks, one prominent at about 
T = 0.5 and a smaller on at about T = 0.1, and see the inset 
for finer detail. The continuous line represents the approx- 
imation furnished by the virial expansion, which obviously 
fails for T < 0.5. 



reader who is mostly interested in the physical insight is 
invited to jump to that subsection. 



B. Specific heat (series expansions) 

The general expression for the pressure (thermal equa- 
tion of state) obtained from the virial theorem is given 
by: 



d(j>ab(rij) 
2d N^' ij dru 



), 



(7) 



where p is the particle number density. The potential is 
a homogenous function of degree — n (Eq. (2)), therefore 



d(j) ab {r) 
dr 



-n(j) ab {r) 



(8) 



Due to this property of the interaction potential we find 
a connection between the pressure and the temperature 
and mean energy: 



P 



T n (U) 



(9) 



This equation is exact for one component and multicom- 
ponent systems and is valid at all temperatures, from 
liquid to solid. 

The next simplification for systems with an inverse 
power inter-particle interaction consists in the depen- 
dence of all excess thermodynamic properties relative to 
the ideal gas on a single density-temperature variable 
[19, 20]. To see why, recall that for a one component 
system the canonical partition function is defined by: 



Zn — 



1 



N\K dN 



exp 



E- 



dr 



(10) 



Here A = /i/(27rmfcsr) 1 / 2 is the thermal de Broglie 
wavelength. The typical distance between particles is 
given by 



l/d 



(11) 



Thus one can use new variables in the integral (10), s* = 
f/l, and the canonical partition function can be rewritten 
as: 



Zn — 



V 



N 



1 



mA dN N N 



exp 



knT 



Q n 



ds 



•N 



where the dimensionless particle number density is de- 
fined by p = ^<J d - This way of writing the partition 
function underlines the existence of the ideal gas contri- 
bution before the configurational integral, and the de- 
pendence of the configurational integral, in the case of 
one component, on a natural parameter, Y = p{-^-^) d l n . 

In the case of a multicomponent system the properties 
of the mixture can be approximated by those of a one 
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component reference fluid [21] with an effective diameter 
defined by: 

a d e =J2^x b a d ab , (13) 

a, b 



where x a = N a /N is the particle number concentration. 
Therefore, the properties of a mixture are defined by the 
effective parameter: 




Nevertheless, in [14] it was shown, that for soft poten- 
tial a more suitable definition of the effective diameter is 
given by: 

o\ = x\g\ + x 2 o\. (15) 

Such a definition leads to a more accurate virial expan- 
sion, as obtained for the present model by [14] using 
molecular dynamics simulations in the temperature range 
0.5<T<5: 



— = 1 + 1.77306r e + 2.362411^ + 2.107981^ + 7.69487T* - 16.23891^ + 27.990871^ - 16.86431^ + 5.469981^ • (16) 

I 



Substitution of the equation (16) to (9) yields the caloric 
equation of state, the specific heat is calculated after that 
by (3). The density corresponding to a given temperature 
is defined as the solution of the equation (16) at P = 13.5. 

Results of the calculations in the frame of the virial 
approach are compared with the simulation results in 
Figs. 3 and 4. Clearly, the virial expansion (16) can- 
not be trusted for temperatures lower than T = 0.5 since 
it was computed from simulations that did not go be- 
low that temperature. This temperature is precisely the 
point at which small micro-clusters become significant, 
and a pure liquid homogenous state description of the 
glassy phase breaks down. Indeed, down to that tem- 
perature the prediction of equation (3) together with the 
virial expansion fits the data excellently well. To un- 
derstand what happens at lower temperatures we must 
await Sect. Ill where the existence of micro-clusters is 
taken explicitly into account. Some comments about the 
limit of temperature going to zero and the relation to the 
Madelung constant can be found in Appendix A 



C. Specific heat (mechanical approach) 

In this subsection we connect the specific heat to the 
bulk modulus of the system. To this aim we begin with 
the microscopic definition of the stress-tensor in an NVT 
ensemble (see, e.g. [22]): 

where pf is the a component of the dimensionless mo- 
mentum of particle i and r°j is the a component of the 
vector joining particles i and j. The first invariant of the 



stress tensor is its trace: 

V(a xx + a yy ) = (£ (Pit + E (P") 2 ) + nU. (18) 

i i 

In order to average (18) one has to take into account that: 

(J2(p7) 2 ) = nm) 2 ) = n-t. (19) 

i 

Thus the average of the first invariant (18) is: 

((a xx ) + (a yy )) = 2pT + np^. (20) 

The pressure is defined as P ~ ((<r xx ) + (<r yy ))/2 and (20) 
yields (9). 

The square of (18) is: 

V 2 ( 

®xx®xx H - ®yyUyy H - ^^xx^yy) (^-0 

=(E (Kf+E (p") 2 ) +2« (e (^) 2 +E (p") 2 ) u 

Wu 2 . 

To compute the average of this equation we need to use 
the fact that ((pf) 2 ) = T and ((pf) 4 ) = 3T 2 . After 
averaging (21) is written as: 

V 2 (((T XX <J XX ) + (VyyVyy) + 2(a xx <r yy )) 

= 4NT 2 + 4N 2 T 2 + ANTn(U) +n 2 (U 2 ). (22) 

The square of (20) is given by: 

V 2 {(a xx ) 2 + (a yy ) 2 +2(a xx )(a yy )) 
= AN 2 T 2 + ANTn(U) + n 2 (U) 2 (23) 
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After substraction of (23) from (22) we have: 



V 



((cr xx a xx ) - (a xx ) 2 ) + {{<J vy a yy ) - (<7 yy ) 2 ) 



+ 2((a xx a yy ) - {a xx ){a yy )) 
= 4NT 2 + n 2 ((U 2 } — {U} 2 ). 



(24) 



This equation has a well defined thermodynamic limit 
since the quantity in square brackets scales like V^ 1 . 
This is seen explicitly in Eq. (25). 

These results allow us now to find exact relationships 
between the specific heat, a caloric quantity, and the elas- 
tic constants which are mechanical quantities. To do so 
we recall the definitions of the elastic constants through 
the stress fluctuations (see, e.g., [23]): 

V 

+ ( C aPiS - Ca/3 7 ,5), (25) 

where C a pjS are the elastic constants and C B p jS are the 
so called Born terms which determine the instantaneous 
clastic constants for any given configuration. 
Substitution of (25) to (24) yields: 



+ V 



\{U 2 ) - (U) 2 ) = 4NT 2 

(C B 4- C B 4- 2C B ) 

\ K - / XXXX ' ^vvvv Ky xxyy ) 



{C XX xx ~t~ C yyyy -\- *2>C xxyy ) 



(26) 



The compression (bulk) modulus in two-dimensional 
systems is: 

K ^ [C X xxx "i" C yyyy + 2C X x yy ) • (^7) 

Recalling Eq. (6) and substituting it and Eq. (27) into 
(26) yields: 



C v K°° - K 

= 1 + 4 

N + n 2 pT " 



(28) 



where the bulk modulus in the infinite frequency limit 
K°° = pT + K B and the Born term is defined by the 
interpaticle interactions [24]: 



This is as much as one can go using exact identities. 
We reiterate that Eq. (28) is very interesting, allowing us 
to connect the bulk modulus to the specific heat. In fact, 
this connection implies that the specific heat measures 
the difference between the bulk modulus and its infinite 
frequency limit. At low temperatures this difference in 
the harmonic approximation as given by: 



K°° ~K = —pT 
4 



(31) 



independent of the solid structure in contrast to the shear 
modulus (cf. [17]). 

The bulk modulus K cannot be computed exactly us- 
ing identities, and we need further information to eval- 
uate it. Fortunately we can estimate the bulk modulus 
from the virial expansion (7) at T > 0.5, since : 



K 



dP 

V 



(32) 



Having done so we can compare the measurements to 
what we expect theoretically. The specific heat as pre- 
dicted from the virial expansion is shown in Fig. 4 as 
the blue (continuous) line. We should stress that com- 
puting Cv from either Eq. (3) or Eq. (28) (using the 
virial expansion (16), yield essentially identical results 
that cannot be distinguished in the blue line in Fig. 4 
for T > 0.5. 



III. SPECIFIC HEAT - THE PHYSICAL 
EXPLANATION 

In this section we propose the physical picture behind 
the existence of the specific heat peaks. We argue that 
the specific heat responds to the micro-melting of the 
clusters - those of small particles at the lowest tempera- 
tures and those of the larger particles at higher temper- 
atures. The large increase in the number of degrees of 
freedom when a particle leaves a crystalline cluster and 
joins the liquid background is the basic reason for the in- 
crease in entropy that is seen as a specific heat peak. We 
can specialize these observations for the model at hand 
(with inverse power potential) or present the discussion in 
greater generality for any model. These two approaches 
are presented in the two following subsections. 



A. Mechanical Equation of State 



k b = — y 



d 2 4>ab{nj) d(f> ab {nj) 



i=£ 3 \ *3 



(29) 



For the present model the "infinite frequency" term (cf. 
[24]) is given by: 



(30) 



In this subsection we employ the mechanical equation 
of state derived above from which the specific heat will 



be computed. To start we define v l u 



and v*i re- 



spectively as the volume of large particle in the whey, 
small particle in the whey, large particle in the solid and 
small particle in the solid. Similarly we denote by e l w , 
e s wl and e s c the energy of a large and small particle in 
the in the whey and in the crystalline phase respectively. 



7 



Needless to say, all these quantities are temperature and 
pressure dependent; we will therefore explicitly use our 
low temperature knowledge concerning vf, and v s c in the 
crystalline phase, but treat the difference v l w — v l c and 
v s w — v s c as constants that we estimate below from our 
simulation knowledge. Similarly we estimate e l c and e s c 
from our knowledge of the hexagonal lattices at T = 0. 
We assume that e l w w and similarly e s w « e s c since 
our simulations indicate a very small change in these pa- 
rameters, see Table I. It should be stressed that the en- 
thalpy change at these pressures are almost all due to the 
PV term. This will result in a semi-quantitative theory 
ascribing the important changes in specific heat to the 
changes in the fraction of particles in curds and whey. In 
other words the number of particles in the whey and the 
number of clusters are all explicit functions of tempera- 
ture and pressure. 



4 


4 


e e 






V'c 


vt 


< 


3.69 


2.07 


3.76 


2.16 


1.43 


0.92 


1.58 


0.94 



TABLE I: Parameters used in the calculation of the specific 
heat 

As the condensed phase consists of clusters of large and 
small partilces, we use the notation N* for the number 
of clusters of n large particles and for the clusters of 
m small particles. In the next subsection we write the 
energy of our system explicitly in terms of these cluster 
numbers. Here however we only need the intensive vari- 
ablcs pi - 2 £„ Ni/N, = 2 £ m K/N p e w = 2Ni/N 
and p s w — 2N^/N which stand for the fraction of large 
particles and small particles in the curds, and large parti- 
cles and small particles in the whey, such that p l c +p^ = 1 
and Pc+ Pw = 1 • Using these variables we can write an 
expression for the volume per particle v = V/N: 

v = + + V -^^P S C ■ (33) 

At this point we need to derive expressions for p c and 
p s c . To do so we need to remember that in the rel- 
evant range of temperatures the large particles in the 
whey can occupy either hexagonal or hcptagonal Voronoi 
cells, whereas small particles can occupy only pentago- 
nal or hexagonal cells [4, 15, 16]. Accordingly there are 
g l w w (2 6 — l)/6 + 2 7 /7 ways to organize the neighbours 
of a large particle in the whey (neglecting the rare large 
particle in heptagonal neighbourhood), but only one way 
in the cluster. Similarly, there are g s w w (2 6 - l)/6 + 2 5 /5 
ways to organize a small particle in the whey. We note 
that this estimate assumes that the relative occurrence 
of the different Voronoi cells is temperature independent. 
While reasonable at higher temperatures [4], at lower 
temperature one should use the full statistical mechan- 
ics as presented in [16] to get more accurate estimates 
of g l w and g s w . This is not our aim here; we aim at a 
physical understanding of the specific heat peaks rather 
than an accurate theory. We thus end up with the simple 



estimates 

pt(P,T) ~ j - g ^ e[{4 _ 4)+P{v e_ vi)]/T > (34) 

Pc(Pi T ) ~ 1 _|_ g« e [(eS-eS,) + P02-<,)]/T ' ( 35 ) 

It is important to note that the combination of Eq. 
(33) together with Eqs. (34) and (35) provides a me- 
chanical equation of state that is alternative to the virial 
expansion presented above. Whereas the latter is best at 
temperature higher than T « 0.5 we expect the present 
one to be best at low temperatures because only the 
present approach takes into account the formation of 
clusters explicitly. The virial expansion by construction 
is a liquid theory. We will now compute C v directly from 
Eq. (28). The peaks in the specific heat are determined 
by the temperature dependence of p c (P,T) and p s c (P,T) 
each which has a temperature and pressure derivatives 
that peaks at a different temperature, denoted as T l {P) 
and T S {P), and see below for details. As said above we 
take Av l = — v l c and Av s = — v s c as approximately 
constants (as a function of temperature and pressure). 
The constants are estimated from the condition that the 
second temperature derivative of p^(P,T) and p s c {P,T) 
should vanish. Explicitly: 

Av e * T e (P*) In g e w /P* , Av s w T S (P*) \ng s w /P* , 

. (36) 

where P* is the pressure for which the peaks in the 
derivatives are observed (13.5 in our simulations). This 
is equivalent to a linear dependence of the specific heat 
peaks as a function of pressure, T e (P)/T e (P*) = P/P* 
and similarly for the small particles. 

In terms of these objects we can rewrite 

v = v c {P,T) + Av e (l- P e c ) + Av s (l-p s c ) , (37) 

To compute the temperature dependence of (§p) T we 
need first to determine its T — > limit, which is deter- 
mined by the first term on the RHS of Eq. (37) as the 
other terms on the RHS decay exponentially fast when 
T — > 0. Since we have already exact results for the bulk 
modulus for the present model, we return to Eqs. (28) 
and (30). We know on the one hand that liniT^o C v = 2 
and that (U)/N w 2.94 over the whole interesting tem- 
prature range, cf. Fig. 3. The compressibility n is re- 
lated to the bulk modulus via k = — (§p) T /v = l/K 
and therefore easily estimated as T — > since there 
(§3?) T w -1/(123.5 - 35T). We use this approximation 
up to T w 0.5. 

Having all the ingredients we can compute C v /N . The 
parameters used were estimated from the numerical sim- 
ulation and are summarized in Table I. Since the aim of 
this subsection is only semi-quantitative, we do not make 
any attempt of parameter fitting, and show the result of 
the calculation in Fig. 5. 
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FIG. 5: Specific heat at constant volume as predicted by the 
simple theory which is based on the mechanical equation of 
state supplied by Eqs. (55), (33) and (34) and (35). Note that 
the theory predicts the two peaks which are associated with 
the micro-melting or micro-freezing of the clusters of large 
and small particles respectively. The magnitude of the peaks 
is too high, reflecting terms missing in the simple approach, 
like the effect of anharmonicity at the lowest temperatures 
which are negative, tending to decrease the height of the low- 
temperature peak. 



Indeed, the theoretical calculation exhibits the exis- 
tence of two, rather than one, specific heat peaks. We 
can now explain the origin of the peaks as resulting from 

the derivatives (jjp^ an d (tUtpJ ■ These derivatives 

change most abruptly when the micro-clusters form (or 
dissolve), each at a specific temperature determined by 
(h^ — hD/lng^ and (h^ — h^/lng^. Note that there can 
be pressures (both upper and lower boundaries) where 
the the sign of (h w — h s c ) or (h l w — h l c ) change sign and 
the peak can be lost. 



B. Caloric Equation of State 

In this subsection we present a more general approach 
which does not take direct input from results derived for 
the inverse power potential. Thus although we use below 
some parameters read from the simulation, the derivation 
is very general any pertains to any distribution of clus- 
ters. To this aim we derive a second equation of state, 
a caloric one. It is quite standard to have two equations 
of state, only for ideal gas and inverse potentials the two 
equations of state degenerate into one. 

Denote the total energy of the system as a sum of E c , 
the energy of the clusters (or curds), and E w , the energy 
of the liquid background (or whey), i.e : 



E = E r 



E,, 



(39) 



transnational, rotational, vibrational and configurational: 

E c = E tr ,c + E rot .c + E v ii,^ c + E con f^ c , (40) 

Ew Ei r w -\- E v *ib w -\- E con f w . (41 ) 

To estimate E tr , c we consider the number of clusters 
of n large particles and of clusters of m small particles 
and write 



(42) 



On the other hand in the whey we follow Eyring [25] and 
Granato [26] and write 



E tr>w =T[Ni + N^]f 



(43) 



where / = 1 — V w s ^ /V w is defined as the fraction of free 
volume in the liquid phase compared to the equivalent 
solid crystalline phase. In other words, 



w i v w u w -r i\ w u w 

V {s) = N e v e + N s v s 



(44) 
(45) 



Similarly, we write 



^rot,c 


n in 


(46) 




= 2T[J2Nin + J2^m] , 

n in 


(47) 




= 2T[N e n + N^](l-f) , 


(48) 


^conj ,c 


= 4J2 N n n + <Y, N ™ m > 

n rn 


(49) 


Econf,w 


= e e N e + e s N s 


(50) 



In terms of these variable we can rewrite 

N n ( n s 

E tr =-T[%- + £- + (pi+p° w )f}, 
z in) (to) 

4 >) (m) J ' 
E vib = NT[2(l-f) + ( P i+p s c )f] , 
N , 



E 



(51) 

(52) 
(53) 



These energies are sums over the degrees of freedom 



Econf = y[4 +e s w + (e* c - 4>c + « ~ Oc](54) 

Summing up all the contributions we need to pay at- 
tention to the order of magnitude of the various terms. 
Since we expect the average size of clusters, at the tem- 
peratures of interest, to be of the order of 30 or so, we can 
neglect safely all the terms that have average size clusters 
in the denominator. With this in mind the expression for 
the energy of the system takes the form 

E o^ro ^ Pc Pc rl I I | s 

N ~ [ 2 /J+ 2 Pc+ 2 Pc 

(55) 
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This is the caloric equation of state that we were after. 
Using it, we can compute C v directly from the thermo- 
dynamic identity 



Cv ^Cp T vu 2 
N N k ' 

where the thermal expansion coefficient is 



a 



i dv 



l r ( w c -vl) ( dpj 
dT 

and the compressibility is 

1 dv \ 
vdPJ T 

U^-O ( dpj 

v l 2 \dP 



(56) 



(57) 



K -v s J ( dp s c 
dT 



+ 



(58) 



{v s c -v s J ( d P s c 
dP 



The last object that we need to obtain for evaluating 

Cy is Cp '. 



Cn 



dE 
df 



P 



dV 
df 



(59) 



Using our expressions (55) and (33) we find 



Cp 
~N 



(2 -Pi 



T 



dp e c 

~df~ 



hi) 



2T 



r 



dp s c 
dT 



(K ~ K) 
2T 



(2-pj-pj) (df 
2 V dT 



(60) 



Having all the ingredients we can sum up the terms in 
Eq. (56). The parameters used were estimated from 
the numerical simulation and are summarized in Table I. 
As before, we did not make any attempt for parameter 
fitting, and show the result of the calculation in Fig. 6. 



C. Discussion 

The bottom line of the simple theory described in 
the previous subsection is that there are two important 
ranges of temperature, first around T * 0.5 where clus- 
ters of large particles begin to form, and a second around 
T « 0.1 where clusters of small particles begin to appear. 
The first important change is also seen in the bulk mod- 
ulus; this is not surprising, since the crystalline clusters 
have a bulk modulus very different from the fluid. Nev- 
ertheless between the clusters we still have appreciable 
fluid regions which act as lubricants for the response to 
shear. We thus expect the shear modulus to change ap- 
preciably only when the small particles begin to cluster, 
in the vicinity of the smaller specific heat peak. We there- 
fore conjecture that the two specific heat peaks are also 
associated with changes in the bulk and shear modulus 
respectively. We expect that any measurements of the 
glass properties connected with bulk and shear moduli 
will show different transition temperatures if these quan- 
tities do not reach simultaneously their K°° counterparts. 

We cannot at this point asses how general is this split 



between bulk and shear moduli, and whether it will be 
seen in generic glasses. We thus leave this point for fur- 
ther research, stressing that we expect this phenomenon 
to appear whenever there exist micro-clusters of preferred 
oredering in the scenario of glass-forming. 

Having an effective equation of state, albeit approxi- 
mate, we can easily compute any thermodynamic deriva- 
tive of interest. We wrote above explicit expressions for 
the compressibility and the thermal expansion coefficient. 
Others are as easily calculated. The point to stress how- 
ever is how non-universal the thermodynamics is. In this 
model we have two specific heat peaks, in others we might 
have one or several. We would also expect a strong pres- 
sure dependence for these peaks. 



IV. FREQUENCY-DEPENDENT SPECIFIC 
HEAT 



By applying time dependent heat fluxes SQ(t) to the 
liquid and measuring the resulting temperature fluctua- 
tions 5T(t), the specific heat can be measured 5Q(t) = 
CST(t). As mentioned in the introduction, measure- 
ments of the specific heat of glassy fluids at low tempera- 
tures can in principle be made under conditions of either 
constant volume (isochoric conditions) or constant pres- 
sure (isobaric conditions), but experimentally isobaric 
conditions are the norm. 

The first and best known measurement of the fre- 
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FIG. 6: Specific heat at constant volume and at constant 
pressure as predicted by the simple theory which is based on 
the caloric equation of state supplied by Eqs. (56)-(59).In 
agreeement with the simulations and the theory based on the 
mechanical equation of state (Cf. Fig. 5) the theory predicts 
the two peaks to both C v and C p which are associated with 
the micro-melting or micro-freezing of the clusters of large 
and small particles respectively. Note that C v > C v as is 
expected from thermodynamics. 



quency dependent complex specific heat was performed 
in glycerol, and we take these experimental results as 
our motivation for this section. We stress from the be- 
ginning that our approach is not particular for glycerol, 
and it can be applied to any other material where, as we 
assume for glycerol, there exists clusters of various sizes 
that determine the dynamical response. In order to de- 
velop a model of the frequency dependent specific heat in 
glycerol we will employ our own model of the glassy phase 
of glycerol. This model assumes that glassy glycerol is 
a heterogeneous fluid on macroscopic timescales. That 
is, that while on very long timescales the liquid phase is 
homogeneous, there exist localized mesoscale domains in 
the fluid that have macroscopic lifetimes. Indeed, inho- 
mogeneities that appear to survive for 10 4 seconds con- 
tribute to the dielectric response in the Fourier domain at 
frequencies as low as 10 -4 Hz in some cases. Clearly such 
imhomogeneities will also contribute to anomalies in the 
frequency dependent specific heat C p (T,lu). We develop 
the theory by deriving expressions for the time-dependent 
enthalpy fluctuations (AH(t)AH(Q)) that are related to 
the frequency dependent specific heat at constant pres- 
sure in terms of the distribution of these heterogeneities. 
The reader is referred to [11] for an introduction to the 
dynamical model of glassy glycerol in which the dielectric 
spectra are computed in great detail. 



A. Frequency Dependent Specific Heat 

By considering a temperature field T(t) = T + 
5T(t),t < 0;T(t) = T, t > and using linear response 



theory on an isobaric ensemble where the appropriate 
Boltzmann distribution is exp(—[3H)/Z, with the en- 
thalpy given by H = E + PV, Nielsen and Dyre [27] 
find that the frequency dependent specific heat is given 
by the form 



C p (7» 



(AH 2 ) 



k B T 2 kWT 2 



B-L Jo 



/ (AH(t)AH(Q))e +iuJt dt. 
Jo 



(61) 

In Eq. (61) AH(t) = H(t) - (H) is an enthalpy fluctua- 
tion away from equilibrium. Therefore we write 



H{t)=Y d N s {t)H s + M l {t)h l 



(62) 



where N s (t) is the number of clusters consisting of s 
molecules in the glassy phase and Mi (t) are the remaining 
molecules in the mobile liquid phase. H s is the enthalpy 
of a cluster of s molecules at a pressure P and tempera- 
ture T 

H S (P,T) = E S (P,T) + PV S (P,T) = (e c +pv c )s + as 2 / 3 . 

(63) 

In Eq. (63) e c (P,T) is the energy /molecule in the con- 
densed phase; v c (P, T) is the volume per molecule in the 
condensed phase; and a(P, T) is the surface energy per 
molecule. Finally hi(P,T) = e/ + Pvi{P,T) is the en- 
thalpy per molecule in the mobile phase. 
Now in equilibrium we can write 



(H)=J2{N s )H B + (M l )h l 

S 

and we also have the sum rule 

J2sN s (t) + M l (t)=M 



(64) 



(65) 



where M is the total number of molecule in the system. 
We can write the enthalpy fluctuations away from equi- 
librium at time t as 

AH(t) = H(t)-(H) = J2(Ns(t)-(N s ))h s = J2 AN ^ 

S S 

(66) 

where 

h s = H s - shi = (e c - ei)s + P(v c - w,)a + as 2 ' z . (67) 

Let us first calculate the equilibrium fluctuations 
((AH) 2 }. To this end we assume that there are no cor- 
relations between the dynamics of clusters of different 
sizes, implying (AN S AN' S ) = 0. Then, using the ex- 
pression Eq. (66) for the enthalpy fluctuations, we can 
immediately write that 



((AH) 2 ) = Y / m 2 )h 2 s 



(68) 



Similarly for the time dependent enthalpy fluctuations 
(AH(t)AH(O)) = ^{AN s (t)AN s (0))h 2 s . (69) 
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We have thus reduced the correlation functions for the 
enthalpy fluctuations into expressions involving the fluc- 
tuations in cluster number for clusters of different sizes 
s. 

To estimate these correlation functions we proceed as 
follows. First we note that the number of molecules N c (t) 
in the clusters can be written as a sum over the clusters 
as N c (t) — J2 S N s (t)s, and consequently the fluctuations 
in the total number of particles away from equilibrium 
are AN c (t) = ^ s AN s (t)s. Then, assuming Gaussian 
fluctuations we estimate the mean square fluctuations of 
the number of particles within some small volume 



<(AiV c ) 2 } a <iV c ) 



(70) 



or rewriting 



]T<(AAgV = £<Ag s . (7i) 



From this equation we therefore see that 
{(AN S ) 2 ) = {N s )/s . 



(72) 



For the time dependent fluctuations therefor, assuming 
an independent Debye relaxation for each cluster, 

(AN s (t)AN s (0)) = (N s )e-V T »/s (73) 

where t s is the lifetime of a cluster of size s. Thus we 
get our final expression for the enthalpy fluctuations in 
equilibrium in terms of the cluster size distribution as 



((AH) 2 ) = Y(N s )(h 2 Js). 



(74) 



and for their time dependent correlations 

W)Aff(O)) = ^(N^e-^ihl/s). (75) 

s 

We now substitute these expressions into Eq. 61 with 
the result 



<w.}«/») 



— l(jJT s 



(76) 



or splitting the specific heat into its real and imaginary 
parts 



(N s )(h 2 Js) 



i x _ (jv s )(^/ s )K) 



kBT 2 ^ 1 + { 



LOT. 



(77) 



We can now use our previous results for the cluster 
distributions in the case of glycerol to find the real and 
imaginary specific heat anomalies in the case of glycerol. 
We do not re-fit any of the parameters used in the calcula- 
tion of the BDS spectra [? ] , we simply use the previous 



Real part of the specific heat times thermal conductivity 




FIG. 7: The theoretical real part of the specific heat 5RC p (lj) 
multiplied by the thermal conductivity for glycerol. 



knowledge at the temperatures indicated, and plot the 
results, fitting only the heat conductivity of glycerol. We 
approximate h 2 s /s ~ (h c — hi) 2 s, such that Eq (76) is 
rewritten as 



C P (T, lo) 



(K - h t f 
k R T 2 



(N a )s 

— IUJT S 



(78) 



Splitting the specific heat into its real and imaginary 
parts 



5RC P (7» 



3C p (T,w) 



(h c - h e ) 2 ^ 

k B T 2 ^ 1 + K) 2 ' 



(h c - hi) 2 



{N s )s{lot s ) 



k B T 2 ^ 1 + (a 



(79) 



(80) 



The resulting curves multiplied by the thermal conduc- 
tivity are shown in Figs. 7 and 8. These should be com- 
pared to Fig. 2 of [3]. The reader can convince himself 
that the theory captures the experimental results quan- 
titatively. 

We would like to stress at this point that the results 
obtained here are equivalent in dynamical contents to 
the computation of the dielectric spectra in [11]. In that 
calculation one focused on the dielectric response e(ui) 
and as here decomposed it into its real and imaginary 
parts e(uj) = 5fte(o;) + i^se(uj). It was found in [11] that 
without the dc contribution (which is absent in the case 
of specific heat) we could write 



»e(w) 



9fe(w) 



E s (Ns)s/[l + (u;r s ) 2 } 
J: s (N s )s(u;t s )/{1 + (ujt s ) 2 } 



(81) 



Once normalized the specific heat spectra are identical 
these spectra . The reason for the identity is in the as- 
sumptions that (m s ■ m s ) ~ s and (hi) ~ s 2 , cf. Eq. 
(67). On the other hand the role of the relaxation times 
t s and the distribution of cluster sizes are exactly the 
same in the two expressions. 
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FIG. 8: The theoretical imaginary part of the specific heat 
QC p (u>) multiplied by the thermal conductivity for glycerol. 



V. SUMMARY AND DISCUSSION 

Probably the most glaring consequence of the calcu- 
lations presented in this paper is that the specific heat 
is a valuable indicator of the interesting physics that oc- 
curs during the glass transition, but this transition is in 
no way universal. The temperature dependence of the 
specific heat is determined by details like inter-particle 
potentials and micro-melting or micro-formation of clus- 
ters. In this sense any hope for universality is untenable. 
Nevertheless we have shown that the specific heat peaks 
herald interesting new physics, leading to fast changes in 
the mechanical moduli which are also associated with fast 
changes in the inhomogeneities that are crucial for the 
glassy behavior, i.e. the formation of micro-clusters. We 
propose that the appearance of two specific heat peaks 
in the case of the binary mixture indicates two different 
ranges for the increase in moduli, the bulk modulus at 
higher temperatures when the first type of clusters form, 
and the shear modulus when the other type of clusters 
form, and the 'lubricating' effect that allows the system 
to shear disappears. All this interesting physics is in- 
dicated by the behavior of the thermodynamics specific 
heat. As for the complex specific heat we have shown, in 
the context of the example of glycerol, that the physics 
revealed by the complex specific heat compared to other 
methods of linear response like Broad Dielectric Spec- 
troscopy are identical. In fact, a straightforward conse- 
quence of our model for glycerol is the prediction that 
the spectra measured from specific heat can be divided 
by the spectra computed, say, from BDS and the result 
should be a constant number. We do not have data for 
exactly the same temperature, but such an experiment 
would be very useful for the near future. 

It is interesting to see in future research whether the 
two specific heat peaks discussed above may be seen in 
other systems, or may be an even richer scenario can 
appear, with more peaks, when more types of clusters 
intervene in the process of glass-formation. 
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APPENDIX A: THE ZERO TEMPERATURE 
LIMIT 

t follows from the simulation results that at lower tem- 
peratures the specific heat is at least close to the value 
of the solid. Therefore, we can write the energy of the 
system in the harmonic approximation: 

dN 

U =U + 22 a iJliQj^ (Al) 

where Uq is the potential energy of the system in the 
reference state, a^- are expansion coefficients and qt is a 
Cartesian coordinate of the deviation of the current posi- 
tion of a particle from equilibrium. There are dN degrees 
of freedom (neglecting translations and rotations of the 
system) in (Al), therefore it follows from the equipar- 
tition theorem that the average potential energy of the 
system in the solid state is [28] : 

<El = ^ + d _ T . (A2 ) 
N N 2 K ' 

Substituting(A2) in (3) immediately yields Cy = 2 for 
two dimensional solids. At the equilibrium configuration 
the potential energy (4) of the binary mixture model (2) 
is given by: 

Due to the scaling properties of the inverse power poten- 
tial it is possible to normalize the interparticle distances 
by the typical distance (11) r. L j — > Sij = rij/l [29] : 

^=cmA (A4) 

where the constant cm — (1/2N) J2ij l/ s "j i s indepen- 
dent of the density. Note that this constant is known as 
the Madelung constant in solid-state physics. Taking into 
account (A4) the average potential energy of a harmonic 
solid (A2) can be rewritten as: 

— = c M p" + -T. (A5) 

The value of the constant cm can be calculated simula- 
tionally using Eq. (A5); the results are shown in Fig. 9. 
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FIG. 9: Evaluated values of the constant cm from the simu- 
lation results. 

One can see that below T = 0.5 the change of this con- 
stant is small. Nevertheless, this change reflects the fact 
that our calculations at the smallest values of the tem- 
perature are not fully relaxed to equilibrium even though 
we took extreme care. Typically at the lowest temper- 
atures the system can be trapped for incredibly long 
times in a local minimum of the energy surface, where 
each local minimum having slightly different equations of 



state [16]. While we expect the Madelung constant to be 
unique for a given crystal, our system here contains clus- 
ters of preferred structures with random orientations [30] , 
and therefore the analog of the Madelung constant is not 
strictly defined. It may very well depend on the protocol 
of cooling. The present best estimate of the value of this 
parameter at the lowest temperatures is cm — 14.649. 

The caloric equation of state (A5) substituted to the 
virial equation (7) gives the following thermal equation: 

^ = (i + f) + 5(c M K)r»/ d . (A6) 

The value of the renormalized constant Cm/cb = 1-394 
can be compared with the result for the two dimensional 
one-component system with hexagonal crystal, which is 
1.268 [31]. The fact that this constant is expected to 
increase in an amorphous solid was anticipated in [29] . 

Finally, we note from (A6) that in contrast to a crys- 
talline solid the thermal (caloric) equation of state here 
remains ambiguous because the value of cm depends on 
the preparation protocol. With this in mind it becomes 
fruitless to seek the anharmonic corrections to equation 
(A6) as in the case of a one component system with a 
well-defined reference state at low temperatures. Never- 
theless we stress that the specific heat at constant volume 
does not suffer from any ambiguity and therefore can be 
taken good as a good indicator of the solidification. 
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